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We consider the interplay of linear double-well-potential (DWP) structures and nonlinear long- 
range interactions of different types, motivated by applications to nonlinear optics and matter waves. 
We find that, while the basic spontaneous-symmetry-breaking (SSB) bifurcation structure in the 
DWP persists in the presence of the long-range interactions, the critical points at which the SSB 
emerges are sensitive to the range of the nonlocal interaction. We quantify the dynamics by devel- 
oping a few-mode approximation corresponding to the DWP structure, and analyze the resulting 
system of ordinary differential equations and its bifurcations in detail. We compare results of this 
£N| ■ analysis with those produced by the full partial differential equation, finding good agreement be- 

tween the two approaches. Effects of the competition between the local self-attraction and nonlocal 
repulsion on the SSB are studied too. A far more complex bifurcation structure involving the 
possibility for not only supercritical but also subcritical bifurcations and even bifurcation loops is 
identified in that case. 
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I. INTRODUCTION 



The studies of Bose-Einstein condensates (BECs) [l|-[3[ and nonlinear optics [J] keep drawing a great deal of 

attention due to experimental advances in versatile realizations of such systems, as well as considerable progress in 

£h ' the analysis of relevant models based on the nonlinear Schrodinger (NLS) -type equations. One of remarkable features 

specific to these fields is the possibility of tailoring particular configurations by dint of suitably designed magnetic 

and/or optical trapping mechanisms (possibly acting in a combination) that confine the atoms in the case of BEC, or 

virtual (photonic) and material structures manipulating the transmission of light in nonlinear optical media. These 

achievements motivate the detailed examination of the existence, stability and dynamical behavior of nonlinear modes 

£> ' in the form of matter or optical waves. The NLS equation [J, Q, as well as its variant known as the Gross-Pitaevskii 

C* (GP) equation [l|-[3| in the BEC context are often at the center of such analysis. 

Within the diverse range of external confinement mechanisms, one that has attracted particular attention is that 
provided by double- well potentials (DWPs). Its prototypical realization in the context of BEC relies on the com- 
bination of a parabolic (harmonic) trap with a periodic potential, which can be created, as an "optical lattice" , by 
the interference of laser beams illuminating the condensate p . The use of a DWP created as a trap for BEC (with 



the intrinsic self- repulsive nonlinearity) has revealed a wealth of new phenomena in recent experiments [7( , including 
the tunneling and Josephson oscillations for small numbers of atoms in the condensate, and macroscopic quantum 
self-trapped states for large atom numbers. Prior to this work, as well as afterwards, motivated by its findings, a 
wide range of theoretical studies investigated such DWP set ting s, including such issues as finite-mode reductions 
and symmetry-breaking bifurcations [a-fla]. quantum effects [16|, and nonlinear variants of the DWP [17| ■ DWP 
; i ■ settings and spontaneous-symmetry-brcaking (SSB) effects in them have also been studied in nonlinear-optical set- 
tings, such as formation of asymmetric states in dual-core fibers |18J |. self-guided laser beams in Kerr media |19j . and 
optically-induced dual-core waveguiding structures in photorefractive crystals [2 Of . 

One of recent developments in both fields of matter and optical waves is the study of effects of long-range nonlinear 
interactions. In the BEC these studies are dealing with condensates formed by magnetically polarized 52 Cr atoms [21J 
(see recent review [22j), dipolar molecules [23(, or atoms in which electric moments are induced by a strong external 



field [24| . Matter- wave solitons supported by the dipole-dipole interactions were predicted in isotropic [25( , anisotropic 
[26|, and discrete J27( two-dimensional (2D) settings, and in the quasi-lD configurations [28|, 1221 (the latter was done 
not only in BEC, but also in a model of the Tonks- Girardeau gas [3(J )■ I n optics, prominent examples of patterns 



supported by long-range effects arc stable vortex rings predicted in media with the thermal nonlocal nonlinearity [31 
as well as the experimental realization of clliptically shaped spatial solitons in these media [32| . 

Our aim in the present work is to examine effects of long-range interactions in the context of DWPs. This is a topic 
of increasing current interest; in the context of dipolar multi-dimensional condensates, it was recently addressed in 
Refs. [331 ] . where the phase diagram of the system was explored, as a function of the strength of the barrier in the 
DWP, number of atoms, and aspect ratio of the system. The possibility of a transition from a symmetric state to an 
asymmetric one, and finally to an unstable higher-dimensional configuration was considered. Here, wc focus on the ID 
setting, and explore different types of long-range interactions, including the dipolc-dipolc interactions, as motivated 



by Rcfs. 28 30j, as well as the interactions with Gaussian and exponential kernels, motivated by nonlinear-optical 
models [34[. In particular, we consider the effect of the range of the interaction, with the objective to consider a 
transition from the contact interactions to progressively longer-range ones. We conclude that the phenomenology of 
the short-range interactions persists, i.e., the earlier discovered SSB phenomena Qj], [lj, [la, [35[ still arise in the present 
context. However, the critical point of the SSB transitions features a definite, monotonically increasing, dependence 
on the interaction range, which is considered in a systematic way. More elaborate scenarios can be detected in the 
case where in addition to the long-range interactions, there is a competing short-range component. In such a case, 
we identify not only the earlier symmetry-breaking bifurcations but also reverse, "symmetry-restoring" bifurcations, 
as well as the potential for symmetry-breaking to arise (for the same parameters) both from the symmetric and from 
the antisymmetric solution branch. Both of these are phenomena that, to the best of our knowledge, have not been 
reported previously. Interestingly, the only example where subcritical bifurcations have been previously discussed in 
the double well setting for GP equations is the very recent one of extremely (and hence somewhat unphysically) high 
nonlinearity exponents in |36J |. 

The presentation is structured as follows. In section II, we present the model and the quasi-analytical two-mode 
approximation, which clearly reveals the system's bifurcation properties. In section III, we corroborate these analytical 
predictions by full numerical results for different types of the long-range kernel. In section IV, we consider a more 
general case, in which the long-range nonlinear repulsion competes with the local attraction. There, we illustrate how 
the competition strongly affects the character of the SSB in the DWP setting. Finally, in section V, summarize our 
findings and present our conclusions. 

II. THE LONG-RANGE MODEL AND THE ANALYTICAL APPROACH 

In the quasi-lD setting, the normalized mean-field wave function ij)(x,t) obeys the scaled GP equation, 



idtip + /J-ip = A^ + s 
where fi is the chemical potential, and 



K (x - x')\il){x')\ 2 dx' 



*!>, (i) 



C = -(l/2)d 2 x + V(x) (2) 

is the usual single-particle energy operator, which includes the confining DWP 

V(x) = (l/2)Cl 2 x 2 + V Q sech 2 (x/W) , (3) 

with il the normalized harmonic-trap's strength; Q <C 1 in a quasi-ld situation in BECs. The nonlinear term with 
coefficient s accounts for the long-range interatomic interactions, s = ±1 corresponding to the repulsion and attraction, 
respectively Note that the contact interaction is not taken into regard in Eq. ([1]), as we aim to focus on the effect 
produced by the long-range nonlinearity (in the gas of 52 Cr atoms, the contact interaction may be readily suppressed 
by means of the Feshbach resonance [2l| ) . In this work, we consider mainly symmetric spatial kernels in Eq. (jTJ) , that 
are positive definite, with the following three natural forms chosen for detailed analysis: the Gaussian, 



the exponential, 



K{x) = -^exp —J , (4) 

K{x) = i-cxp (-M) , (5) 

and the cut-off (CO) (alias generalized Lorentzian) kernel, 

K(x) = —o- 3 (x 2 +o- 2 r 3 / 2 . (6) 

The width of the kernels, a, determines the degree of the nonlocality All three kinds of the kernels go over into the 
(5-function as a approaches zero, in which case Eq. (TTJ) turns into the usual local NLS/GP equation. All the kernels are 
normalized as the 5-function, i.e., J_ K(x)dx = 1, and the norm of the wave function will be used in the usual form, 
N = J_ \i/j(x,t)\ 2 dx. In what follows below, we adopt typical physically relevant values of the scaled parameters, 

namely, Ct = 0.1, Vq — 1 and W = 0.5, in which case the two lowest eigenvalues of linear operator C with potential 
([3]) are numerically found to be luq — 0.1328 and ui\ = 0.1557. 



A. The two-mode approximation 

The spectrum of the underlying linear Schrodinger equation (s = 0) consists of the ground state, with wave function 
uo(x), and excited states, ui(x) (I > 1). In the weakly nonlinear regime, the Galerkin-type two-mode approximation is 
employed to decompose the wave function tp(x, t) over the minimum basis constituted by the ground state Uq and the 
first excited state ui, associated to the eigenvalues ujq and Wi, respectively. For this purpose, it is more convenient to 
use a transformed orthonormal basis composed by wave functions centered at the left and right wells, viz., {ip^, tp R } 
= {(uq — ui)/v2, (ua + Ui)/v2}. Without the loss of generality, ipl,r are both chosen to be positive definite. Thus, 
the two-mode approximation for the wave function is defined as 

ip(x, t) = c L (t)ip L (x) + c R (t)ip R (x), (7) 

where Cl and c R are complex time-dependent amplitudes. Substituting this into Eq. ([1]), we obtain 

ic L ip R + ic R (p R = (ttc L - [ic L - u>c R )ip L + (ttc R - pc R - u>c L )ip R 

+ s\c L \ 2 (cnp L + c R tp R ) / K(x - x')ip 2 L (x')dx' + s\c R \ 2 (c L ip L + c R tp R ) / K(x - x')(p R (x')dx' 



(8) 
+ s [( c 1 c r + \ c L\ 2 c R )ip L + (c* L c R + c L \c R \ 2 )(p R ] / K(x - x')ip L (x')f R (x')dx', 



where the asterisk and overdot stand for the complex conjugate and time derivative, while Q = (loq + u>\)/2 and 
u = (cjj — ujo)/2 are linear combinations of the two lowest eigenvalues. Next, we project Eq. ((5]) onto the single-well 
states lpl.fi.1 which involves the following overlap integrals: 



Vo =11 K{x - x')<p 2 L (x')tp 2 L (x) dx'dx, 
?/i = / / K(x — x')ip 2 L (x')ip 2 R (x) dx'dx, 



(9) 
m = 1 1 K(x - x')ip 2 L (x')(f L (x)ip R (x) dx'dx, 



m =11 K i x - x')tp L (x')(p R (x')(p L (x)ip R (x) dx'dx. 



The equations hold if subscripts L and R are swapped, or the variable x and x' are interchanged, due to the symmetry 
of kernel K. In Fig. Q] we show the values of the four integrals, 770,1,2,3, as functions of parameter a, for the three 
types of kernels defined in Eqs. Q-©. We note that 772 and 773 remain negligible for any value of cr, and 770,1 arc 
both positive when the kernel K is positive definite. Naturally, when a is small, which corresponds to the limit of 
a nearly local nonlinearity, 770 is much larger that the other three overlap integrals, due to the weak overlapping of 
the single-well states ^pl,r- As a increases, 770 decreases while 771 increases and, finally, they tend to become equal. 
Regarding the values of these integrals, the upcoming analysis is conducted in two situations: 770 much larger than 
all others, and 770 being on the same order of magnitude as 771. Some value of the ab is to be fixed as a boundary 
between the two cases. Since 770 is always large, the key point is to set up a rule for comparing 771 with 772,3. To this 
end, we define 77 rc i = 7/1 — maxd^l, I773Q and the criterion is stated as follows: if ?7 rc i > 0.01, 771 is taken into regard 
in the analysis of the two-mode approximation; otherwise, 771 is insignificant, and only 770 is kept. For the kernels 
considered in this context, the criterion yields <Jb = 2.96, 1.56 and 1.91 for the Gaussian, exponential and Lorentzian 
kernels, respectively. 

Thus, with a < a^, all integrals 771,2,3 are omitted, and only 770 is retained. Then, the projection of Eq. ((SJ) onto 
the two-mode set, ifiL,R, leads to the following ODE system: 

ic L = (0 - n)c L - uc R + si] \c L \ 2 c L , 
ic R = (Cl- fi)c R - uc L + st] Q \c r \ 2 c r . 

We then introduce the Madelung representation, cl, R = Pl,r£ ' R with real time-dependent pl,r and #l,_r, and 
derive from Eq. (jlOp a set of equations for pl,r. and 9l,r- 

Pl = up R sin 9 

a ( m , P R a 2 (H) 

9l = (ii — \i) + uj — cosw — ST] p L 
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FIG. 1: (Color online) The overlap integrals 7701 ?7i, V2 an d 773 as functions of width a of the kernels of the following types: 
Gaussian (left), exponential (middle), and Lorentzian (right), as defined in Eqs. Q-©. 

where 9 = 9j, — 9r is the relative phase between the two modes, and the equations for pr and Br are obtained by 
interchanging subscripts L and R, and 9 with — 9, in Eq. (jTTj) directly. We focus on steady solutions to this system, 
i.e. pl.r = 9l,r = 0. Then, for solutions with nonzero amplitudes, 9 may only take values or tt, which correspond, 
respectively, to equal or opposite signs of real stationary solutions for cl and Cr. Through a straightforward algebra, 
three stationary solutions are thus found: the symmetric solution, with 9 = and p\ R = (p,— u)o)/srjo, existing when 
p > lo q (p < cj ) for s = 1 (s = —1), i.e., for the repulsive (attractive) long-range interactions; the antisymmetric 
solution, with 9 = ir and p\ R = (p — uii)/sr)o, existing when p > ui\ (p < us\) for s = 1 (s = — 1); and an asymmetric 

solution, with p\ R = (s(p — Cl) ± y/(p — Cl) 2 — 4uj 2 )/2rjo- As we assume 770 > 0, for the repulsive case (s = 1) the 
asymmetric state exists only if 9 = n, i.e., it bifurcates from the antisymmetric solution when p > 51 + 2cj. On the 
contrary, in the attractive case, it emerges from the symmetric state (9 = 0) when p < Q — 2a; . These conclusions 
agree with the general principles of the SSB theory, according to which the attractive/repulsive nonlinearity breaks 
the symmetry/anti-symmetry of solutions with equal numbers of particles in the two wells [12|, [lj, Ua,l35j. 

Next, we consider the other situation, in which both rjo and r\\ are taken into consideration, while the other two 
arc neglected (i.e., a > er&). In this case, projecting Eq. ([HJ onto (fL,R results in the following system: 



ic L = (fl - p)c L - ujcr + sc L (\c L \ 2 r) Q + |c_r| 2 t7i), 
ic R = (O - p)c R - uc L + scr{\cr\ 2 t]o + |c L | 2 ?7i). 

In terms of the Madelung representation, we transform Eqs. (fT2|) into 



(12) 



p L = upr sin 9 



PR 



9 L = (p-n)+to—cos9 

PL 



svopl ~ smP% 



(13) 



with the equations for pR and 9r produced by swapping subscripts L and R, and 9 with —0, as before. In this case, 
a set of three stationary states are again obtained: the symmetric one, 9 = and p L R = (p — u>o)/s(r)o + 771) ; the 
antisymmetric state, 9 = tt and p\ R = (p — Ux)/s{t]q + -qi). Finally, the asymmetric state, 9 = n (9 = 0) for s = 1 
(s = — 1) and p\ R = ((/i — fl) J sr] ± ^J \p — Cl) 2 /rj 2 . — 4o; 2 /At7 2 )/2, where A77 = 770 — 7/1 > for the kernels we consider, 
exists at p > Cl + 2uj7]o/Arj (p < Cl — 2uj7]o/Ai]) for s = 1 (s = —1). The value of p at which the SSB bifurcation 
happens is tagged as p cl . Since rjo/Arj is an increasing function of a, the value of p CI increases (decreases) as a grows 
for s = 1 (s = — 1). 

The two-mode approximation is a powerful means for identifying different steady states of the underlying problem. 
The solutions found above are also used as initial conditions in solving the full NLS equation ([1]), as reported in 
section IIIII 

B. The bifurcation analysis 

The two-mode approximation strongly facilitates the qualitative analysis of the SSB bifurcation, as well as exploring 
the system's dynamics. To proceed, we define the population imbalance between the two wells, 



z = (N L - N R )/N = (\c L \ 2 - \c R \ 2 )/N, 



(14) 



where Ni ji = \cl,r\ 2 = p 2 L R , hence the total norm is TV = Nl + TVr. Together with the relative phase, 6 = 6l — 6r, 
we eventually derive the following dynamical equations: 



i = 2cj v 1 — z 2 sin 9 

■ 2uzcos9 Ar ( 15 ) 

= -== — srjNz. 

V 1 — z 2 

This form of the equations is relevant for both cases, when only 770 or both 770 and r\\ dominate, as discussed before. 
Note that 77 stands for 770 in the former case, and for A77 in the latter one. Equations (|15[) take the Hamiltonian form, 

~m < 16 » 



d 



Z 



with Hamiltonian 



H = 2u\J\ - z 2 cos 6 - -si]Nz 2 . (17) 

Equations (|15[) possess the stationary solutions, (zi,#i) and (22,^2), with z\ = Z2 = 0, B\ = 0, #2 = 7T, which 
represent the symmetric and the antisymmetric solutions, respectively. Besides that, the asymmetric stationary 
solutions may exist when TV > TV cr , with TV cr = \2uj/t]\, taking the form of 

Auj 2 

e = ir(e = o), z 2 = i-—, (is) 

for s = 1 (s = —1). The asymmetric solution emerges from the antisymmetric (symmetric) one through a pitchfork 
bifurcation, in the case of the defocusing (focusing) nonlinearity. Since 77 = 770 for a < erf, and rj = 770 — 7/1 for a > erf,, 
Fig. [1] suggests that 77 is, generally, a decreasing function of a, and consequently 7V cr is increasing with respect to a. 
This way, the bifurcation takes place at larger value of TV when the the nonlocality range is wider, which is consistent 
with the conclusion concerning /u cr . 

Next, from Eq. (|15[) we derive the equation of motion, 



z = -Wz + \ij\NzyjAu 2 - icu 2 z 2 - z 2 , (19) 

which leads to the system 



z = p, 



(20) 



\p= -4w 2 z + \ti\Nz^Auj 2 - Auj 2 z 2 - p 2 . 
When TV < TV cr , there is a unique stationary solution, p = z = 0, which is a fixed point of the center type. When 



TV > TV cr . the origin, (0, 0), becomes a saddle, with another pair of fixed points (centers), p = 0, z = ±4 1 — -r—r^ — 5-, 

V \m n 2 

representing the asymmetric solutions. Figure [2] shows the phase space of system (|20| . along with the linearization 
near the fixed points for the example of the Gaussian kernel with a = 5.0 and TV = 0.5, in which case TV cr = 0.33. 

III. NUMERICAL APPROACH 

We first consider the repulsive interaction case (s = 1). Branches of stationary solutions to the full partial differential 
equation ([1]) are explored for all three kernels defined in Eqs. (|4])-([6]) with various values of parameter a . The results 
are plotted in Figs. [3]-|H where the solutions arc expressed in terms of TV, the number of atoms, as a function of 
the chemical potential fi. The stationary solutions are sought by employing a fixed-point Ncwton-Raphson iteration 
onto a finite difference scheme with Ax = 0.1 and using the continuation of the solutions with respect to fi. The 
linear stability of each solution i))q, is analyzed by considering the standard linearization around it in the form 
ip(x,t) = V'o + e[a(x)e xt + b*(x)e x ']. The relevant linear eigenvalue problem is written as 

VaHQ- (21) 
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FIG. 2: (Color online) Phase diagrams in the Gaussian-kernel model with with a — 5.0 and N = 0.5, in which case N CT — 0.33. 
The left panel is obtained from complete system (|20[) . displaying three fixed points, (0,0) and (±0.7498,0). The middle and 
right panels are the phase planes of the linearized system near fixed points (0,0) (middle) and (0.7498,0) (right). 



where operators L\, Li are defined as 

1 r°° 

-dl + V-v + s / K {x - x')\i)o(x')\ 2 dx' 

K (x - x')ipo(x')^{x')<j){x') dx' 



Litf> = 

L 2 4> = s 



K(x- x')il>o(x')il>o(x f )(l>(x f ) dx', 

-OO 



(22) 



for any function <j>. The stationary state is called unstable if there exist any eigenvalues A with 5R(A) ^ 0, otherwise 
it is stable (i.e. all corresponding eigenvalues are purely imaginary). 

We present the results of Gaussian kernel in detail as an example. In each panel (Fig. |3|), the solid blue line with 
highest value of N for any fj, among the three branches is the symmetric stationary solution; the continuation of it 
to the linear limit (N —> 0) shows that it starts from fx = u>q (the eigenvalue associated to the ground mode of the 
underlying linear system), and it is stable for any fj,. The branch with a part of it denoted by dashed red line is the 
antisymmetric solution, arising from the linear limit at it, = uj\, corresponding to the first excited state. It starts as 
a stable state from the linear mode (denoted by solid blue line) and, as /x increases, after some critical point /j cr , 
it is destabilized due to the emergence of the asymmetric branch through a supercritical pitchfork; notice that the 
asymmetric branch remains stable after it occurs. This bifurcation was predicted by the two-mode analysis in the 
previous section. 

The four panels of Fig. [3] are the bifurcation diagrams of the model with the Gaussian kernel, for a = 0.1, 1, 5 and 
10. The main parts of the panels display the numerically found solutions of the full system, while the small plot in 
each panel shows the corresponding analytical solutions predicted by the two-mode approximation, as obtained in 
Section [TTJ The numerical and analytical branches of the solutions demonstrate good agreement in all the four cases. 
Recall that we have obtained the critical values at which the bifurcation takes place, fi cr (or equivalently N cr ), in an 
analytical form. As for the Gaussian kernel (with the chosen border at <7b = 2.96), a = 0.1 and a = 1 are categorized 
as belonging to the first case, with solely tjq taken into account, among all the overlap integrals. The analytically 
predicted value is /i cr = 0.1671 (since /i cr = ft + 2w), while its numerically found counterpart is /j cr = 0.1684 for both 
a = 0.1 and 1. On the other hand, a = 5 and a = 10 pertain to the the second case, in which both 770 and r\\ are kept. 
In this case, the analytical prediction is /u cr = Q + 2ujr)o/ A77, which depends on 770,1 and thus varies for different values 
of a. The predicted values of /z cr is 0.1757 and 0.2120, for a = 5 and 10, respectively, while the respective numerical 
values are /j cr = 0.1745 and 0.2075. Despite a small discrepancy between the two groups of the values (numerical 
versus analytical), the analytical results still succeed in predicting the trend of the behavior of /i cr , viz., /u cr increases 
with the growth of a, or, in other words, the SSB bifurcation takes place at higher values of N cr , as shown in Fig. |U 

Next we consider cases with the long-range interaction based on the other two kernels, viz., the exponential and 
Lorentzian ones, given by Eq. (|S|) and ([B]), respectively The corresponding bifurcation diagrams, similar to the case 
of the Gaussian kernel, are presented in Fig. [5] and Fig. |6l with the same notations as adopted above. 

Lastly, we also briefly discuss the self-attractive case, with s = — 1 in Eq. |T]). Two examples are selected, using 
the Gaussian kernel with a = 1 and 5, to help illustrating the similarities and differences with the self-repulsive case. 
The complete bifurcation diagrams are plotted, for this case, in Fig. with the notation similar to that introduced 
above. 

In this case, the two branches arise from the linear modes at /i = u>o and u±, pertaining to the symmetric and 
antisymmetric stationary solutions, which exist for /_j < Wo or fj, < oj\, respectively. This time, the supercritical SSB 
pitchfork bifurcation occurs on the symmetric branch, leading to the emergence of the asymmetric state. The pitchfork 
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FIG. 3: (Color online) The normalized norm, N, of the numerically found stationary solutions of the underlying GP equation 
with the Gaussian kernel, for the case of repulsive interaction (s = 1), as a function of chemical potential \x. The four panels 
correspond to cases with different values of a, viz., a = 0.1 (top left), a — 1 (top right), a = 5 (bottom left), and, finally, a = 10 
(bottom right). The stationary solutions predicted by the two- mode approximation for each case are shown in the small plot in 
the top left corner of each panel. The blue solid lines and red dashed lines denote stable and unstable solutions, respectively. 
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FIG. 4: (Color online) The critical value of the chemical potential, fi CT , at which the supercritical pitchfork bifurcation takes 
place in the case of the self-repulsive nonlinearity with the Gaussian kernel, as a function of a. The blue solid line and the red 
dashed-dotted line denote the numerically obtained values and their counterparts predicted by the two-mode approximation, 
respectively. Note that there is a small jump at a — 2.96 on the dashed-dotted line, due to the definition if the border between 
the two situations [described by Eqs. (|10[) and (|12[) . respectively] in the two-mode analysis. 



is found at /i cr = 0.1213 and 0.1128 for a = 1 and 5 respectively, in good agreement with the critical values predicted 
by the two-mode approximation, which are 0.1212 and 0.1104. In each panel of Fig. [3 the diagram produced by the 
two-mode approximation, corresponding to its numerical counterpart, is displayed in the top right corner, exhibiting 
a very good agreement between the two. 
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FIG. 5: (Figure online) Norm N of the numerically found solutions and their counterparts predicted by the two-mode approx- 
imation (shown in the corners) for the case of the self-repulsive nonlinearity with the exponential kernel, as a function of /i. 
Here and in the next figure, the value of a is chosen as a — 1 (left), a — 5 (right), and the notation is the same as in Fig. [3] 
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FIG. 6: (Color online) The norm of the numerical and approximate analytical solutions for the case of the self-repulsive 
nonlinearity with the Lorentzian kernel, as a function of /i. 



IV. COMPETING SHORT- AND LONG-RANGE INTERACTIONS 



In this section, wc consider a model illustrating effects of the competition of the long-range interactions with local 
ones, based on the following GP equation, 



id t ip + jiip = I/ip + s 



K(x-x')\ip(x')\ 2 dx' 



V^ + slVf^, 



(23) 



where coefficient g accounts for the contact nonlinearity. The interactions are fully attractive or repulsive if s and g 
are both negative or positive, respectively. However, in such cases the results turn out to be similar to those reported 
above for the nonlocal nonlinearity. A more interesting case arises for sg < 0, i.e., for the competing interactions. 
Effects of the competition on ID solitons were recently studied in detail in Ref. |29| . 

We use the analysis presented in the previous section as a starting point for our considerations here. Thus, we fix 
the long-range interaction coefficient s = ±1, and then vary the local interaction coefficient, g from to — 1 (g < 
implies the local attraction), using Gaussian kernel (j4]) with a = 5 and s = 1 as a case example. As the attractive 
local interactions grow stronger, interesting phenomena emerge in the development of the corresponding bifurcation 
diagram, as is shown in Fig. [8] 

When g is close to zero, the three branches of stationary states are similar to those in the case of pure long-range 
interactions, as shown in the bottom left panel of Fig. [3J An example for g = —0.2 is presented in the top left panel 
of Fig. [51 where the norm N of each branch is larger at any value of /x, in comparison to the case of g = 0. Also, the 
bifurcation takes place at a larger critical value, /x cr = 0.1805. 

As g grows more negative, we notice that the symmetric and the antisymmetric branches, while monotonically 
increasing with /i at small values of N , switch to become monotonically decreasing functions of fi at a large value 
of N. Simultaneously, we observe the bifurcation of two asymmetric states, each from a different branch. While it 
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FIG. 7: (Color online) Norm N of the numerical solutions, together with their analytical counterparts (shown in the small plots 
on the top right corners), for the case of the self-attractive nonlinearity (s = —1), as a function of /i. The Gaussian kernel is 
taken, with a — 1 (left) and a = 5 (right). The notation is the same as in Fig. [3] 
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FIG. 8: (Color online) Norm N of the numerically found solutions in the case of the competing interactions for Gaussian kernel 
with a = 5, as a function of fi. The coefficient of the long-range interactions is s — 1, while the local-interaction coefficient 
is g — —0.2 (top left), —0.3 (top right), —0.35 (bottom left) and —0.45 (bottom right). The notation is the same as in Fig. [3] 



is commonly known that attractive and repulsive interactions favor bifurcations from symmetric and anti-symmetric 
branches, respectively, here we encounter the first (to our knowledge) example that features bifurcations from both 
branches. This unusual situation extends up to g = —0.32. The complete bifurcation diagram is displayed in the top 
right panel of Fig. [FJ for the case of g = —0.3. In this case, the symmetric and antisymmetric branches arise, as usual, 
from their linear limits at fj, = wo and n = lo\ , respectively. A supercritical pitchfork bifurcation occurring (as before) 
on the antisymmetric branch gives rise to an asymmetric state at (i = —0.195, destabilizing the antisymmetric one. 
What is, however, different here is that this asymmetric state "survives" only within a narrow parametric interval 
before it merges back into the antisymmetric branch through another subcritical ("symmetry-restoring") pitchfork 
at fi — 0.3263. This bifurcation loop is reminiscent of that reported previously for two-component solitons in local 
ID and 2D models with competing self- focusing cubic and self-dcfocusing quintic nonlinearities [37| . For higher /z, 
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another asymmetric branch arises from the symmetric one through a subcritical pitchfork at (i = 0.3465; it is especially 
important to highlight that this pitchfork, which destabilizes the symmetric branch, is a subcritical one, hence the 
asymmetric state is unstable too. It is observed that, after its emergence, the norm of the asymmetric branch decreases 
as /i decreases (dN/dn > 0), before starting to rise at the fold point at \i = 0.2768, featuring dN/d/j, < after that. 
The stability of the state changes at the turning point, which is in agreement with the well-known Vakhitov-Kolokolov 
(VK) criterion [38| , according to which the slope of the branch determines its stability. In Fig. [9j a blowup provides 
a clearer view of these bifurcations. 
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FIG. 9: (Color online) The left panel: a segment of the top right panel in Fig. [8] for g — —0.35, clarifies details of the 
bifurcation picture. The right panel: the asymmetry-measuring intergal characteristic, J_ °° x\ip\ 2 dx, for the antisymmetric 
solution and the asymmetric one arising from it, as a function of fi, in the same case as in the left panel. The loop illustrating 
the asymmetric solution emerging from (through the symmetry-breaking supercritical pitchfork) and merging back (through 
the symmetry-restoring subcritical pitchfork) into the antisymmetric branch. The notation is the same as in Fig. [3] 



The symmetric and the antisymmetric branches continue to increase their norms with the growth of /z, each of them 
featuring its own turning point, at n = 0.391 and /i = 0.3853 respectively, and change their monotonicity thereafter. 
However, this change on the antisymmetric branch is not accompanied by a stability change. To demonstrate that 
this is not an exception to the VK criterion, we define the linearization operators, 



L— = L\ — L-2 



■3 ff |V| 2 , 



(24) 



using notation £1.2 defined in Eq. (|22|). It is known from Ref. |39j that the instability arises, when the slope condition 
of the VK criterion is violated, if \n(L + ) — n{L_)\ = 1, n standing for the number of negative eigenvalues of each 
operator. The state under consideration is gencrically unstable if \n{L + ) — n{LJ)\ > 1. In our case, the antisymmetric 
branch has n(L-) = 1 (due to its single-zero-crossing configuration profile, which is a zero mode of £-). Meanwhile, 
n(L+) = 1 is true before the branch turns to the left, and consequently \n(L+) — n(L-)\ = 0. Hence, even though 
the slope condition is violated (dN/d/j, > 0), the relevant theorem suggests that the VK criterion does not actually 
apply here (and no stability change occurs). After the turning point, the stability is expected since n(L+) = 2, and 
the slope condition is satisfied, which is in accordance with our observations. 

Next, as the local interactions continue to grow stronger, another regime is explored, as shown in the bottom left 
panel of Fig. [51 with g = —0.35 as an example. In this case, the phenomenology of the symmetric branch is similar 
to that in the previous situation; however, the two sequential bifurcations taking place on the antisymmetric branch 
no longer arise. Hence, the entire antisymmetric branch remains dynamically stable. The asymmetric branch still 
emerges from the symmetric one through a subcritical pitchfork, and later becomes stable past the fold point due to 
the change of slope, in accordance with the VK criterion. 

The above regime persists until g = —0.38. After that, the pitchfork bifurcation occurring on the symmetric branch 
switches from subcritical to supercritical, leading to the emergence of a stable asymmetric branch, as shown in the last 
panel of Fig. [5] This shape of the bifurcation diagram, consisting of the three branches, persists as far as g decreases 
to —1. Note that during this process, the turning points of the symmetric and antisymmetric branches keep getting 
lower (with respect to N), finally leading to a situation where the turning points do not exist, and the two branches 
immediately go left starting from their linear limits, i.e., the symmetric and the antisymmetric solutions only exist 
for fi < ujq, lui when g < —0.69, —0.6, respectively. 



11 

V. CONCLUSIONS AND FUTURE CHALLENGES 

In this work, we have presented a systematic study of the interplay between the long-range nonlinear interactions (of 
cither sign, repulsive or attractive) and linear double-well potential (DWP) in the ID setting. The two-mode approx- 
imation has been developed, that accounts for the nonlocality but retains the structure similar to that established 
before in the case of the contact interaction. This conclusion demonstrates that the fundamental phenomenology 
of the spontaneous symmetry breaking (SSB) in the DWP persists in the presence of the longer-range interactions, 
although the critical properties themselves (e.g., the critical values of the chemical potential and norm), at which 
the SSB bifurcation occurs, giving rise to the asymmetric steady states, are sensitive to the precise interaction range. 
In particular, our analysis has revealed a monotonous increase of the critical values as a function of the interaction 
range, in the case of the self-repulsion. A considerably more elaborate phenomenology was revealed by the analysis in 
the context of competing long- and short-range interactions. In that case, a delicate interplay between the strengths 
of the nonlocal repulsion and local attraction gives rise, in addition to predominantly attractive and predominantly 
repulsive regimes, to mixed ones with complex bifurcation phenomena. On the one hand, symmetry-breaking effects 
were shown to arise from each of the relevant solution branches; in some cases, they are accompanied by reverse 
bifurcations, to form closed loops. On the other hand, some of the previously studied supercritical bifurcations, such 
as the one occurring on the symmetric branch, could become subcritical, being subsequently coupled to additional 
fold bifurcations. All of these effects are absent, to our knowledge, in models with a single cubic nonlinear term, being 
consequences of the competition. 

These results may be a motivation for studies in a number of future directions. In particular, it would be relevant to 
consider the generalizations of the DWP setting to higher dimensions, such as, the four-well configuration, which was 
recently demonstrated to yield a much richer phenomenology in the case of the contact interactions |40f . In the same 
connection, it is relevant to note that the few-mode approach, used in this work for the consideration of the existence 
and stability of the symmetric, antisymmetric, and asymmetric states, as functions of the interaction strength, may 



be applied to other structures, such as multi-dimensional bright [25| . dark [4l[ and vortex [3l| solitons, especially in 
cases where such structures are stabilized by nonlocal nonlincaritics. A natural objective of such an analysis may be 
to unravel effects of the interaction range on structural stability of nonlinear waveforms. Such studies are presently 
in progress and will be reported elsewhere. 
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